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Abstract 



We perform a numerical investigation of the thermodynamics during the collapse of 

a charged (complex) scalar field to a Reissner-Nordstrom (RN) black hole in isotropic 

coordinates. Numerical work on gravitational collapse in isotropic coordinates has recently 

shown that the negative of the total Lagrangian approaches the Helmholtz free energy 

F = E — TS of a Schwarzschild black hole at late times of the collapse (where E is the 

f»o ' black hole mass, T the temperature and 5* the entropy). The relevant thermodynamic 

CNJ • potential for the RN black hole is the Gibbs free energy G = E — TS — &hQ where Q is 

~4^. ' the charge and $# the electrostatic potential at the outer horizon. In charged collapse, 

there is a large outgoing matter wave which prevents the exterior from settling quickly to a 

CN ' static state. However, the interior region is not affected significantly by the wave. We find 

numerically that the interior contribution to the Gibbs free energy is entirely gravitational 

and accumulates in a thin shell just inside the horizon. The entropy is gravitational in 

origin and one observes dynamically that it resides on the horizon. We also compare the 

numerical value of the interior Lagrangian to the expected analytical value of the interior 

Gibbs free energy for different initial states and we find that they agree to within 10—13%. 

The two values are approaching each other so that their difference decreases with more 

evolution time. 
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1 Introduction 

Recent numerical studies of gravitational collapse in isotropic coordinates, both in 4 + 1 dimen- 
sions (5D) [TJ [2] and 4D [2], have provided numerical evidence that the negative of the total 
Lagrangian at late times of the collapse process approaches the Helmholtz free energy E—TS 
of a Schwarzschild black hole, where E, T and S are the mass, temperature and entropy of 
the black hole. For a stationary black hole, T = Kk,/2tt and S = A/4h where k is the surface 
gravity at the horizon and A is the horizon area. Though T and S both contain h, their 
product TS does not and this opens the possibility for a classical investigation of black hole 
thermodynamics via the free energy [1]. There is in fact an argument based on the Euclidean 
action that relates the negative of the total Lagrangian to the free energy of a stationary black 
hole [H [3] . This association can be tested numerically during gravitational collapse by tracking 
the Lagrangian and comparing its numerical value at late times to the expected free energy 
from standard black hole thermodynamics. This was carried out numerically for the first time 
for the collapse of a 5D Yang-Mills instanton [TJ to a Schwarzschild black hole in isotropic 
coordinates. A numerical study in isotropic coordinates was then carried out for the collapse 
of a 4D and 5D massless scalar field to a Schwarzschild black hole [2]. These works constituted 
a classical numerical study of black hole thermodynamics; quantum mechanics does not enter 
the picture and Hawking radiation is not observed. Prior studies of thermodynamics during 
gravitational collapse have consisted mainly of analytical or semi-analytical work on black hole 
entropy and Hawking radiation during shell or dust collapse [4]-|ll| and the approaches have 
been quantum mechanical or semi-classical. 

For the Reissner-Nordstrom (RN) black hole, the relevant thermodynamic potential is the 
Gibbs free energy G = E — TS — ^hQ where &h is the electrostatic potential at the horizon 
and Q is the charge of the black hole O [12]. In this paper, we investigate numerically in 
isotropic coordinates the thermodynamics during the collapse of a charged scalar field to a 
RN black hole. Among other things, we study the association between the total Lagrangian 
and the Gibbs free energy. Charged collapse leads to a large outgoing matter wave and the 
exterior takes a long time to settle to a static state. However, the interior spacetime (inside the 
outer horizon) is hardly affected by this wave and settles more quickly than the exterior. One 
can show analytically that the interior Gibbs free energy Gi n t is equal to —TS. The interior 
contribution of the negative of the total Lagrangian can be obtained numerically and then 
compared to Gj n j. At late times of the collapse, we find that the matter Lagrangian tends 
towards zero in the interior; the interior Lagrangian stems entirely from the gravitational 
sector. It is negative and accumulates just inside the horizon where the spacetime is not 
static. In short, the entropy of the charged black hole is gravitational in origin and stems from 
the dynamical interior near the horizon. That the entropy stems from the interior is in accord 
with analytical work on the canonical quantization of the RN black hole [10] where the entropy 
was obtained via explicit counting of microstates in the dynamical interior (see [UJ for the 



uncharged case) . Our numerical results for the negative of the gravitational Lagrangian in the 
interior agrees with the analytical value of Gi n t to within 10—13% depending on the profile of 
the initial state. There are sharp changes in the gradients of the metric and matter functions in 
the near-horizon region and this places limits on how far in time one can evolve before the usual 
monitors of the code, such as the ADM mass, begin to deviate from their conserved values. 
The numerical graph representing the Lagrangian and the analytical graph representing the 
free energy approach each other with time and it is clear that given more time evolution the 
difference between them would continue to decrease. Of course, to increase the time evolution 
one needs to increase the resolution and one reaches a limit where the computing time is no 
longer practical. As a consistency check, we also implement a procedure for prolonging the 
evolution time in the exterior region. The outgoing matter wave disperses and the metric in 
the exterior is observed to approach the RN static exterior metric. Unlike Schwarzschild, the 
matter Lagrangian now makes a contribution to the exterior due to the presence of a static 
electric field. 

The coordinate time t in isotropic coordinates coincides with the time measured by an asymp- 
totic observer at rest. The viewpoint of the static asymptotic observer is appropriate for 
studying black hole thermodynamics. The temperature of a stationary black hole is the tem- 
perature as seen by an observer at rest at spatial infinity |13j and its entropy represents a 
measure of an external observer's ignorance of the internal configurations hidden behind the 
event horizon |14[ [T5] . In particular, the association between the total Lagrangian and the free 
energy of a black hole that is found in numerical studies of collapse in isotropic coordinates will 
in general not hold in other coordinate systems. Unlike the action, the Lagrangian depends on 
the choice of time coordinate i.e. on the foliation of the spacetime. For example, one should 
not expect the Lagrangian in numerical studies of gravitational collapse in Painleve-Gullstrand 
(PG) coordinates [16] to be associated with the free energy of the black hole. The reason is that 
in PG coordinates the coordinate time represents the proper time of freely-falling observers 
[17j not asymptotic observers. This is also true of numerical studies of charged collapse that 
have been carried out in "Eddington-Finkelstein" [18] and double-null coordinates [191 [20] as 
there is again no coordinate representing the proper time of an asymptotic observer in these 
cases. 

Our paper proceeds as follows. We first express the exterior RN metric in isotropic coordinates. 
We then derive the equations of motion in those coordinates: the wave equation governing 
the complex (charged) scalar field, Maxwell's equations for the gauge field and Einstein's 
equations governing the metric field. We then discuss how the initial states are obtained using 
the shooting method. Integral and differential expressions in isotropic coordinates are then 
derived for the conserved charge Q and mass M. During the evolution these quantities should 
remain constant and this allows one to monitor the accuracy of the code at each time step. 
Expressions for the gravitational and matter Lagrangian as well as the interior Gibbs free 
energy are then derived. We finally present the thermodynamic results from the numerical 



simulation. 

We work in geometric units G = c = 1 where energy, mass and time have units of length. 
Coulomb's constant is set to unity so that electric charge has units of length also. The metric 
has signature (— 1, 1, 1, 1). Spacetime indices are in Greek and run from to 3. 

1.1 Reissner-Nordstrom metric in isotropic coordinates 

During our numerical simulation in isotropic coordinates, we expect the metric to settle to the 
Reissner-Nordstrom (RN) metric at late times. We therefore need to express the RN metric in 
isotropic coordinates. A general spherically symmetric time-dependent AD metric in isotropic 
coordinates takes the form [214 [22], [23] : 

ds 2 = -N(r, tfdt 2 + ip(r, tf{dr 2 + r 2 dtt 2 ) (1) 

where ip(r, t) is referred to as the conformal factor and N(r, t) is called the lapse function. Note 
that the isotropic radial coordinate r is not the areal radius. We assume asymptotic flatness 
so that N and ip are both unity at infinity. Note that the coordinate time t coincides with the 
time measured by a clock at rest at infinity. Our goal is to find the analytical expressions for 
iV and ip that correspond to the RN metric. In standard coordinates, the RN metric is given 

by HZ! 

, 9 / 2M Q 2 \ l9 / 2M Q 2 \~ l , /9 n ,„9 

ds 2 = -(l- — + ^)dt 2 + (l-— + ^) dr' 2 +r' 2 dn 2 (2) 

where M and Q are the mass and charge of the black hole and r' is the areal radius. The 
function f(r') = 1 - 2M/r' + Q 2 /r' 2 has zeroes at r ± = M ± \/ M 2 — Q 2 where r' + and r'_ are 
referred to as the outer and inner horizon respectively. Matching the metric ([1]) to the metric 
© yields ip 4 = r' 2 /r 2 , N 2 = f(r') and the equation ^ = jpj$^. With the condition that 

i/j and N are unity asymptotically, the latter has solution r' = r + M + (M 2 — Q 2 )/Ar so that 
the exterior region of the RN metric in isotropic coordinates is given by 
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The outer horizon in isotropic coordinates is situated at r + = \J M 2 — Q 2 /2. At this location, 
the lapse function is zero: N(r+) = 0. The metric ([3]) does not cover the interior region of 
the RN black hole; it covers the static exterior region twice. The minimum value of r' = 
r + M + (M 2 - Q 2 )/Ar is r' + = M + sj M 2 - Q 2 ; in metric © both the region r > r + 
and r < r + correspond to the exterior region r' > r' + . The Killing vectors in the interior 



region between the two horizons of the RN black hole are all spacelike [2] and the spacetime 
is nonstationary in this region. To cover the interior region of the RN black hole in isotropic 
coordinates the functions N and if) in (TjQ) must be time-dependent. In standard coordinates, 
the spacetime becomes nonstationary when one crosses the outer horizon into the interior 
because the function f(r') switches sign and the radial coordinate becomes timelike and the 
time coordinate becomes spacelike. In isotropic coordinates, the metric coefficients in (TjQ) are 
positive-definite; they do not switch sign but instead become time-dependent in the interior. 
During our numerical simulation, the metric functions N(r, t) and ^(r, t) in the isotropic metric 
(TjQ) are nonstationary in the region r<r + , reflecting the true nature of the interior spacetime. 
It is only in the exterior region r>r + that the metric approaches the static form ([3|) at late 
times. 

We do not observe the second (inner) horizon or timelike singularity associated with the RN 
metric ([2]); we observe one horizon at r + (where N{r + ) = 0) and a spacelike singularity as 
in the Schwarzschild case. This is in accord with the findings of previous numerical work on 
charged collapse |18l [TUI |2"U] . The inner horizon of the RN metric ([5]) is an artifact of exact 
staticity (and exact spherical symmetry) [T7] and it has been known since pioneering work in 
the 90's [24], that the inner (Cauchy) horizon is unstable to perturbations. 



2 Evolution and constraint equations in isotropic coordinates 

2.1 Matter sector 

For matter, we consider a complex (charged) scalar field coupled to an electromagnetic field 
Au,. The matter Lagrangian density C m has a local U{\) gauge symmetry and is given by [25] 



Cm = ~\ (x-,v + »eA M x) <T (X;u - ieA u x) - -^i^ F^ (4) 

where a semi-colon denotes covariant differentiation evaluated with metric ([1]), a bar denotes 
complex conjugation and F^ = A U] ^ — A^ v is the electromagnetic field tensor. Spherical 
symmetry reduces the number of gauge components from four to two: only At = Aq and 
A r = A\ are non-zero. Gauge freedom allows one to further eliminate A r . This leaves At as 
the only non-zero component and for simplicity we label it a. The matter fields are therefore 
X = x(r, £) and a = a(r,t). Lagrange's equations of motion for matter are 

V«^ 7T— = (5) 

dq- a dq 

where q is a generic field. 



2.1.1 Equations of motion for scalar field \ 

Applying ([5]) to the scalar field x yields the wave equation 

X-, »v 9^ + ieAp gT (2xj u + ieA„x) + ieA K „ gTx = • 
With spherical symmetry, the three terms in the above equation reduce to 
1 
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where x — dx/dt and %' = dx/dr. Equation (jfij) now reads 
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d r {Ni; 2 r 2 x>) - ^ (X + *eo X ) = . (7) 
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For numerical purposes, we split the above second order equation into two first order evolution 
equations. To this end we define the quantity 



P = jjr (x + ieax) 



In terms of p, equation (j2|) simplifies to 



p = ~2 d r {Nip r x') ~ ieap. 



Rearranging (|5J) we obtain 



N 
ft 



X = lfi\P- 



ieaip 6 x 

N 



Equations (j9]) and (fTOl) are the evolution equations for p and x respectively. 



(8) 

(9) 
(10) 



2.1.2 Equations of motion for gauge field a 

Applying ([5]) to the gauge field A^ yields Maxwell's equations: 
1 



2tt 



F mp 9 up + iex {X;im ~ ieA^x) ~ i&X (x-,u + ieA^x) = 0. 



(11) 



Spherical symmetry reduces the above to two equations, one for fi = t and fj, = r. The equation 
for fj, = t is 

a" N'a! 2i\)'a! 2a' 2irieN ,_ = 

where p was defined in ([8]) . The above can be rewritten as 



4 + i^4-^r-;^: + -^(^-xp) = o (12) 



1 fip 2 a'r 2 



-Ji dr ( — ^r~ ) = 27rie (xp - xp) ■ (13) 



The above is a constraint equation for the gauge field a (one can view it as Poisson's equation). 
Maxwell's equation for fi = r yields 

a! Naf 2ipa' _, , 

which can be expressed as 

g = 2-KieN^ 2 (xx - xx) (15) 

with 

9=—. (16) 

Equation (|15p is the evolution equation for g. The equation governing the gauge field a is not 
an evolution equation but an ordinary differential equation: a' = g N/ip 2 . 

2.2 Gravitational sector 

2.2.1 Stress-energy tensor 

The stress-energy tensor T^ v appearing in Einstein's field equations can be calculated from 
the matter Lagrangian @ . We first express the latter in terms of the new quantities p and g 
defined in ([8]) and (fTB]) respectively: 

c _w__i^_ + J_ (17) 

Lm ~ 2^ 2 2V> 4 8vn/> 8 ' [ ' 

The energy-momentum tensor is given by 

T at 3 = ~ 2 ^ZaJ + C ™ 9a/3 (18) 



and its non-zero components are 



T ee = „ ( -7s ~ x' x' + i — 7T ) ?" 2 2"W = T 90 sin 2 (9 
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tt 2 VV> 12 V> 4 4vrV 6 

2.2.2 Field equations 

Einstein's field equations are given by G^ = k 2 T^ where k 2 = 8tt G = 8ir and G^ is the 
Einstein tensor evaluated with metric (P). The G rr equation yields evolution equations for 
the conformal factor tp while Gu and G r t yield constraint equations. G$e yields an ordinary 
differential equation for the lapse function N (G^ yields the same equation). 

The G rr equation is 

pr = (2ri{j' 2 N 3 - Arip 2 N^ 4 + 2^'N 3 ^ + 2rNipip 5 - 2r4>N^ + 2rN'^'N 2 ^ 
rw I N 6 V 

2 /.,„ „2 X (20) 



+N ' w V) = y(f + ^-4^) 



We split the above second order equation into two first order evolution equations. For this 
purpose we define 

* = -£■ < 21 > 

K is in fact the negative of the trace of the extrinsic curvature for spacelike hypersurfaces at 
constant time t. Equation (I20p can now be expressed as 



K K 2 6^' M 1\ 3N' (2i\)' 1~ 
TV " T ~ IF \J + r) ~ NiA \~J~ + r, 

3/£ (XX_ , PP ff 



4 V V' 4 V' 12 4tt^ 8 
Equations (|21|) and (|22|) are evolution equations for ip and K respectively. The Gu equation is 

4 fo/2/3 n,/A r 2 Ar2 ,/A k 2 N 2 f pp x'x' , ff 2 



and using K can be cast in the form 



^ NW) = |(^ + |1 + _|L)_^. (23) 



Equation (|23p is the energy constraint. The G r t equation is 



(fal/N - ^ f N^ + ^JVty) = ^ 2 — ^ (XP + X'p) 



^ 2 iv v r ; 2^ 

which can be expressed as 

K' K 2 
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K' K 2 

= + x'p) ■ (24) 



3 4-0 
Equation (j24|> is the momentum constraint. The G00 equation gives 

Urip 5 ipN - 8r^ij 2 N - 4rip 5 ^N + N'ip 2 N 2 - 2rip' 2 N 3 + 2rt/;'N 3 tp 



^ 2 N 3 

+2riVW + r*W) = y (f " x¥ + ^ ) ^ 



Eliminating the time derivatives using the G rr equation ([20]) . the above equation reduces to 
an ordinary differential equation for N 



i/j 2 r \ril) 3 ) N r \rxlA) " V 4 4?r0 



A ^ + „* 1^)=— 5^ + x^ (25) 



The evolution equation (|22p for X can be made more numerically stable by combining it with 
(1231) to obtain 



K _ 2 6^ /^ 1\ 3JV' /2^ 1 

= (2^ + n/>") V- 

n/> 51 ; 8vr^ 8 



(26) 



To summarize, the evolution equations for the matter sector are (jlOp and (JSj) for the pair (x,p) 
and ()15p for g respectively. The evolution equations for the gravitational sector are (|2ip and 
P6P for the pair (ip, K). The gauge field a and lapse function N obey the ordinary differential 
equations (|16p and (|25p respectively. There is one constraint equation in the matter sector, 
namely "Poisson's" equation f)13[) . In the gravitational sector, the energy and momentum con- 
straint are given by equations (|23p and (|24j) respectively. Once the initial states and boundary 
conditions are specified, the evolution of the fields is unique. 



3 Initial states 

Let xi an d X2 be the real and imaginary part of the complex scalar field \. We begin by 
choosing the initial configuration for Xlj X2, Xi an d X2- ^ is convenient to pick a configuration 
such that the momentum constraint (|24|) is trivially satisfied i.e. where the right hand side of 
p4p is zero so that K' = initially. Since K = asymptotically, this implies that K is initially 
zero everywhere. We use three different profiles that satisfy this property: 

I: Xi=X2 = j0^ Xl = ~X2 = ^t 

II: Xi = X2 = A 2 (tanh (Ai - r) + 1) X 'i = ~X2 = ^f 

III: Xl =V^-z ; X2 = xi=X2 = 0. 

(Af+r 2 ) 

Profile I// represents a real scalar field with no charge. One is free to choose the values of 
the different Aj and these in turn determine the values of the conserved mass M and charge 
Q. The initial states for the metric functions tp and iV and for the gauge field a are obtained 
by solving three coupled second order differential equations: the energy constraint (f23|) . the 
ordinary differential equation ([25]) and "Poisson's" equation fjl3[> . These can be split into six 
first order differential equations: 

y/W (x'x! (X + ieax) (x - ieax) C 2 \ 




N 2 4vrr 4 ^ 8 y 

k 2 (x + ieax) {x ~ ieax) k 2 



3x¥ 




AN 2 4ip 4 



r 2 ip 2 
where the last three equations introduce the new variables A, B and C: 

The functions range from to a large distance R, the outer computational boundary repre- 
senting "infinity" . The spacetime is asymptotically flat and boundary conditions at r = R that 



are consistent with this are N(R) = 1, ip(R) = 1 and B(R) = 0. Gauge invariance allows us to 
set a(R) = 0. The initial configuration is non-singular at r = and this implies that A(0) = 
and C(0) = 0. To solve equations (|27p with the above boundary conditions, we use a shooting 
method with a point located between the origin and R. From the boundary conditions, two 
fields are known at the origin and four fields are known at R. We make an educated guess 
as to the value of each field at the missing end (a total of six numbers.) We then evolve the 
fields from both ends toward the middle point using fourth order Runge-Kutta. The goal is to 
modify the six unknown numbers until the fields match in the middle. This can be achieved 
in a rather straightforward way using a six dimensional Newton method |26j . The system of 
linear algebraic equations can be solved using an lower/upper matrix decomposition |26] , A 
first estimate for the fields can be obtained by solving the equations for x = which applies to 
a scalar field. The fields can be made to fit at the middle point to within a few parts in 10 15 . 
We then increase x by a small amount and use the preceding values of the fields to estimate 
the new ones. We repeat the procedure until a high enough charge to mass ratio has been 
obtained. It is worthwhile to note that if the fields match at the middle point, it can be shown 
that their derivatives also match. 



4 Expressions for the conserved charge Q tot and mass M tot 
in isotropic coordinates 

In this section we obtain expressions in isotropic coordinates for the total (conserved) charge 
Qtot and mass Mtot- During the simulation, these quantities should remain constant and this is 
used to monitor the simulation. The mass M and charge Q of the black hole itself is significantly 
less (magnitude wise) than Mtot and Qtot respectively because a considerable amount of charge 
and mass is expelled during the collapse process. M and Q will be evaluated in the results 
section with a different expression. For the analytical stationary RN metric ([2]) or Q, there is 
no outgoing matter wave so that M and Q are equal to Mtot and Qtot respectively. We leave 
them expressed in terms of M and Q because in the results section they are evaluated with our 
black hole values of M and Q (not with our values of M to t or Q to t as these do not correspond 
to the mass and charge of the black hole respectively) . 

The total conserved charge Qtot is given by the general formula |17j 

Qtot = ^j F a ?dS aP (29) 

where S is a closed two-surface bounding the charge distribution and dS a g is the two dimen- 
sional surface element given by [T7] 

dS a p = -2n [a r p] ^d 2 6 . (30) 

10 



n a and r a are the timelike and spacelike normals to the two-surface and the square brackets 

denote antisymmetrization: n\ a rffi = {n a rp — npr a )/2. The induced surface element on S 

is yfad 2 9 where a= det(a ab) where (jab is the induced metric on S (A,B is either 8 or </>). 

For the metric ([I]), we obtain n a = (— A", 0, 0, 0) and rp = (0,i/j 2 ,0,0). The only non-zero 

components of F a ^ are F tr = —F rt = jjQtkx - The only non-zero components of dS a p are 

dStr = —dS r t = N ip e r 2 sin 9d6d<p. Putting these results together we obtain an expression for 

the total charge: 

V> 2 rV , , 2 | 

Qtot - -^^\ r=R -a H \ r=R (31) 

where we used the fact that asymptotically, in the large R (infinite) limit, N and tp approach 
unity. Using the constraint equation fjl3[> . one can express the charge Qtot in the integral form 

f R 

Qtot= 2irie (xp- XP) r 2 dr. (32) 

Jo 
From the above expression, one can define an effective charge density 

P= ~2 (XP-XP) ■ (33) 

The total (conserved) ADM mass Mtot for an asymptotically flat spacetime is given by [T7] 

1 /",. . . ^ j2/ 



Mtot = -—<t(k- k )^d 2 6 (34) 

Sn Js 

where S is the two-sphere at spatial infinity, k is the trace of the extrinsic curvature of S 
embedded in St, the three-dimensional spacelike hypersurface at constant t. k$ is the trace of 
the extrinsic curvature of S embedded in flat spacetime. k and ko can be readily evaluated 
and one obtains k — ko = 4?// /ip 3 . The mass M to t is equal to 

M to t = -2rV Vi =iJ = -2R 2 ip'\ r=R (35) 



where we used that ip = 1 asymptotically. Using (|23p one can express the mass in the integral 
form 

M k 2 f R f x 'x' PP g 2 2K 2 \ 2 5 

Mtot = Tj \^ + W 2 + ^¥~^) r ^ dr - (36) 

5 Free energy 

There are arguments, based on the Euclidean action [3], for relating the negative of the total 
Lagrangian of a stationary black hole to its free energy. As far as we know, to date, there is no 

11 



direct analytical proof of this; it has only been tested numerically during gravitational collapse 
to a Schwarzschild black hole [IJ [2] . One of the goals of this paper is to test this association 
numerically for the case of charged collapse where the relevant thermodynamic potential is the 
Gibbs free energy [21 [12]. To do so, we first need to determine an expression for the Lagrangian. 

5.1 Gravitational and Matter Lagrangian 

The gravitational action Sg [g] is a sum of the Einstein-Hilbert term Sjj [g] , a boundary term 
Sb [g] and a nondynamical term Sq [T7] : 

Sg [g] = S H [g] + S B [g] + S = J ' L G dt (37) 

where Lq is the gravitational Lagrangian and 

S H [g} = -^Jv^gRd 3 xdt. (38) 

The Ricci scalar in isotropic coordinates is given by 

R = — Tpnr (I8rip 2 Nib 3 + GribNib 4 - 8tb'N 3 - 4rt//'iV 3 (39) 

rN d ip b 

- 6rNipip 4 - rN"N 2 i) - 2rN'^'N 2 - 2N'N 2 ^) . 

The Ricci scalar is not the Lagrangian density as it contains second derivatives of the metric 
tensor. Besides the Hilbert action Sn[g], one also requires the boundary term Ss[g] to obtain 
Einstein's field equations [T7j- So ensures that the action Sg is finite and equal to zero in flat 
spacetime. We do not need to spell out Sb (or S* ) to determine the gravitational Lagrangian. 
As in pQ, Lg can be obtained simply by integrating out by parts the second derivative terms 
appearing in R and this yields 

L G = - ! R (srViVV + 8r 2 i\y 2 - ^f^ ) dr . ( 4 Q) 



4 7 V N 

Note that Lq contains only first derivatives of the metric functions. The matter Lagrangian 
is given by the integral of (fT7l) : 

7 3 A f R ( pp x'x' . g 2 



L m = l£ m ^-gd*x = 4n J (£j,-^ + ±jj)^\N\r 2 dr. (41) 



/o V 2 ^ 12 2^ 4 8vr^ 
The total Lagrangian is the sum L m + Lg- 
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5.2 Gibbs free energy of the Reissnner-Nordstrom black hole 

We can easily calculate the Gibbs free energy G = E — TS — &hQ of the RN black hole (E 
is the ADM mass M, T and S are the temperature and entropy respectively and &h the 
electrostatic potential at the horizon). The temperature is given by T = Hk/2tt where k is 
the surface gravity of the black hole and S = A/ Ah where A is the area of the event horizon. 
For a RN black hole, k = \/M 2 — Q 2 /r' + and the area of the horizon is A = 4-7rr^_ , where 
r' + = M + \/M 2 — Q 2 is the (areal) radius of the outer horizon. Therefore the product T S 
is equal to y M 2 — Q 2 /2. The electrostatic potential at the horizon is given by &h = Q/ r '+ 
[3j [12] . The Gibbs free energy of the RN black hole is equal to 

G = M- l -^M 2 -Q 2 Q = 

2 V M+ v / M 2 -Q 2 (42) 

= \VM 2 -Q 2 . 

For Q = 0, expression (J42]) reduces to the (Helmholtz) free energy M/2 of a Schwarzschild black 
hole [U[2]. Note that though h appears in T and S, it does not appear in the product TS and 
hence does not appear in the expression (J42]) for the Gibbs free energy. 



5.2.1 Analytical formulas for the exterior and interior contribution to the free 
energy 

There is an interior and exterior contribution to the gravitational Lagrangian Lq, the integral 
given by (|40p . The analytical expressions for ip and N in the exterior static region can be 
extracted from metric ([3]): 



M 2 -Q 2 



N = i an i • ju 



4r 
M M' 1 -Q 



M M 2 - Q 2 



1/2 

(43) 



Inserting the above functions into (|40p and integrating from an isotropic radial coordinate of 
zero to the outer horizon r + = \J M 2 — Q 2 /2 yields the exterior contribution Lc ext to the 
gravitational Lagrangian 

i roc . . i 

L Gext = j (8rV^V + 8r 2 A^' 2 ) dr = --{M + ^M 2 - Q 2 ) . (44) 

There is no charge residing in the exterior of the RN black hole; there is an electric field but 
no scalar field. The contribution to the matter Lagrangian (]4ip in the exterior stems entirely 
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from the electromagnetic part, 

1 f'°° n 2 O 2 

Lm -L ^-Nr 2 dr = % . (45) 

M ^ 2j r+ V 2 2(M + ^M 2 - Q 2 ) V ' 

The exterior contribution to the Gibbs free energy is the negative of the total exterior La- 
grangian, 

G ext = -L Gext - L Mext = VM 2 - Q 2 . (46) 

The interior contribution is therefore equal to 

1 



G int = G- G ext = -- VM 2 - Q 2 = -TS = -r+ (47) 



where G is the Gibbs free energy of the RN black hole given by (|42h . As in the Schwarzschild 
black hole [U [2], the interior contribution to the free energy is equal to the negative of the 
product of temperature and entropy i.e. Gi n t = —TS. We will see that the negative dip in 
the interior occurs in a thin shell near the radius of the horizon, as in the Schwarzschild case 
[23]. Our numerical results in the interior will be compared to the analytical expression (|47p . 



6 Numerical Results 

6.1 Code and initial state 

We work in geometrized units where G = c = l. We also set Coulomb's constant to unity. In 
these units radius, mass, energy, time and charge have units of length. The coupling constant 
e has units of inverse length. Its value is arbitrary (there is no experimental value for it 
in our model). We set it equal to unity and leave the scale unspecified (one can work in 
centimetres, metres, etc.). One unit of time is the time it takes light to travel a radial distance 
of one unit in flat spacetimqj. For the spatial grid, we introduce a new coordinate x such that 
r = 2a;/ (1 — x). We use a constant grid in x, with x ranging between and 1. This maintains a 
larger concentration of points in the interior than the exterior. The spatial and time increment 
are taken to be Aa; = At = 5 x 10 -5 . The fields that are governed by an evolution equation are 
evolved via a fourth order Adams-Bashforth-Moulton (ABM) scheme [26]. The fields N and 
a are governed by an ordinary differential equation and are obtained by iterating backwards 
starting from infinity using equations ()25[) and (|13p and their boundary conditions respectively. 
This is performed at each time step using again ABM. The simulation was carried out for the 
three different initial profiles discussed in the section on the initial states. The results for 
profiles I, II and III are presented in tables 1, 2 and 3 respectively. M and Q represent the 



1 If the scale were specified to be 3 x 10 5 km, then one unit of time would correspond to one second. 
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mass and charge of the black hole itself, whereas Mf t and Qtot represent the total (conserved) 
mass and charge. We will see shortly how M and Q are evaluated. — Lj n t and Gj n j represent 
the interior Lagrangian and Gibbs free energy at late times respectively. The discussion that 
follows will focus on one initial state with profile II. This will allow us to highlight the essential 
features of the thermodynamics that are common to all cases. The initial state for profile 77 
amounts to specifying the values of three constants Ai, A2 and A3. We will discuss the middle 
case in table 2: Ai = 2, A2 = 0.11 and A3 = 0.032. The total (conserved) ADM mass and 
charge corresponding to this initial state is M to tai = 1-241 and Qtotai = —0.842 respectivel)o 
Note that this differs considerably from the mass M = 0.721 and charge Q = —0.215 of the 
black hole itself because a significant amount of mass and charge are expelled in the collapse 
process. The total mass and charge are useful for monitoring the accuracy of the code and are 
checked at each time step. The ADM mass begins to deviate from its original value before the 
total charge. For the above initial state, the numerical results and plots are obtained up to 
an evolution time of i = 23.5, just below 5% ADM mass deviation. The ADM mass is known 
to be very sensitive to tiny deviations in the metric and matter functions [2] i.e. the functions 
can be evolving well even when the ADM mass begins to deviate from its original value. For 
example, the numerical curve for the metric function ip continues to approach the analytical 
curve of the RN exterior metric even at 5% deviation from the ADM mass. 



6.2 Plots of metric and matter functions, mass and charge of black hole 

We made plots of all the relevant functions. The evolution of the conformal factor ip and lapse 
iV are shown in figures [2] and [3] respectively. We denote r + (t) as the isotropic radius where the 
lapse A^ crosses zero at a time t. Black hole formation is taken to coincide with A^ crossing zero 
[211 I22t [T| since radially outgoing (or ingoing) null geodesies in the interior (r < r + (t)) do not 
cross the 7V = spacelike two-surface as it evolves (see figure Q]). At very late times the N = 
spacelike two-surface is identified as the apparent horizon of the stationary RN black-hole (for 
more details sec [23]). The radius r + (i) increases with time and has a numerical value of 0.344 
at t = 23.5; this is the radius of the apparent horizon. In the interior at late times, the lapse N 
shows almost no dependence on r and approaches zero asymptotically from below (i.e. N — > 
as TV — > 0). The conformal factor ip has its peak value at r + . In the interior, it also has no 
spatial dependence and approaches zero asymptotically from above i.e. ip — > as ip — > (see 
figure E|). As in [18} [j~9l [20] we do not observe an inner horizon; the spacetime structure in 
the interior resembles that of the uncharged case [2] . The absolute value of the complex scalar 
field x is shown in figured! It approaches a constant value inside r + and the only region where 
the time or radial derivatives are high is near the radius r + . Note the outgoing matter wave. 



2 Switching the sign of xi and X2 in the initial state leads to a positive charge. The sign has no bearing on 
the thermodynamics. 
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In the RN solution, the electromagnetic-field tensor is given by F tr = Q/r' 2 [17] where r' is 
the areal radius. A short calculation yields g(r) = Q/r 2 where g is given by (|16p and r is 
the isotropic radius. This, in combination with the relation r + 2 = (M 2 — Q 2 )/A, allows to 
determine the charge and mass of the black hole: 

+ 1_ (48) 

M = ^/Ar + 2 + Q 2 . 

The above are evaluated at the latest time. 

6.3 Electrostatic potential and charged shell 

The field a is shown in figure Inside r + , it is almost perfectly constant in space and its 
value decreases with time. Between r + and the position of the outgoing wave, it decreases 
smoothly as r increases. On the other side of the outgoing wave, it decreases faster and this 
is consistent with the fact that the outgoing wave contains charge. This is what would be 
expected from basic classical electromagnetism for a situation where two charged shells are 
concentric and the outside one is moving outward. It can be seen from figure [6] that the 
charge density of the black hole becomes increasingly concentrated near r + as time progresses. 
At late times it is concentrated just inside the apparent horizon (located at isotropic radius 
r = 0.344). In short, at late times one has a charged shell [23] with a constant potential a 
inside. The value of —a extracted at the apparent horizon in our numerical simulation at late 
times is —0.174. The value of &h = QA+ is —0.153 where r' + = M + \/M 2 — Q 2 is the areal 
radius at the horizon. The percentage difference between —a and &h is around 12%, which is 
explained by the fact that the outgoing wave hasn't reached infinity yet. Note that &h should 
be compared to —a since $# = —A^^\h = —M\h [32] where \h means evaluated at the 
horizon. Here £** is the Killing vector of the stationary RN black hole and in our coordinate 
system it is given by £** = (1, 0, 0, 0). Recall that a = At in our numerical simulation. What is 
important here is that the charge accumulates in a shell at the radial location of the horizon 
|23] and has an electrostatic potential which is numerically close to &h- This provides us 
with a clear physical interpretation of the work term 5W = §h$Q in the first law of black hole 
thermodynamics: it is the work needed to bring a charge 5Q to the charged shell located at the 
horizon. Note that though the charged shell is located at the radius of the horizon, the proper 
radial distance between the shell and the origin decreases towards zero with time because the 
conformal factor ip approaches zero in the interior. The charged shell is effectively "collapsing 
towards the center" (see [23] for details). 
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6.4 Lagrangian and free energy 

6.4.1 The interior region 

Our goal here is to compare the negative of the total interior Lagrangian —L (gravity + matter) 
to the interior Gibbs free energy G{ n t- The first thing to note is that the matter Lagrangian, 
plotted in figure [9l makes basically zero contribution to the interior at late times. The in- 
terior Lagrangian is entirely gravitational. The integral of the negative of the gravitational 
Lagrangian from the origin to a given r is shown in figure [71 There is a clear negative dip just 
inside the apparent horizon that becomes thinner and closer to the horizon with time. Since 
Gint = —TS this suggests that the entropy of the charged black hole accumulates at the horizon 
and that it is gravitational in origin. In the interior, the metric is not static near the horizon 
(see figure ED but dynamical so that the entropy is associated with the dynamical interior in 
accord with some previous analytical studies |11[ PT0| [To] . This was also observed numerically 
for the Schwarzschild case [HE]- The numerical value of the minimum in figure [7J is plotted in 
figure [10] as a function of time together with the interior Gibbs free energy Gi n t = —r + . Note 
how the two graphs approach each other with time. At the late time t = 23.5, the gravitational 
minimum (which we label — Lj n j) is —0.306 and Gj n j = —0.344, for a difference of 10.9% (see 
table 2). It is clear from the graphs that the two values would continue to approach each other 
if the code could evolve further in time. The ADM mass has deviated by 5% at this point. 
This is due to sharp changes in the gradients in the near-horizon region. 

We verified that the grid size determines how far the simulation runs before the ADM mass 
begins to deviate from its original value. This is shown in table 4. We used the following 
sequence for the number of points: {1250,2500,5000,10000,20000} with the corresponding 
step sizes shown in table 4. The times at which the ADM mass deviated by 1% and the 
times by which it deviated by 5% are listed. The table shows clearly that decreasing the step 
size increases the time evolution. The numerical results and plots presented above are for 
the highest resolution Ax = 5.0 x 10 -5 and are taken at t = 23.5 just below 5% ADM mass 
deviation. Going beyond this resolution increases significantly the computing time so that it 
no longer becomes practical. 

The results for all three profiles with different initial states are shown in tables 1, 2 and 3. Note 
that it is possible to have the total charge \Qtot\ be greater than the total mass M tot (see table 
1). However, in all cases, \Q\ is less than the mass M and there are no naked singularities. The 
discrepancy between — Lj ni and Gj n j is shown for all cases. As one can see, the discrepancy 
depends slightly on the initial state. However, in all probability, this difference is not due to 
fundamental physical reasons. It is simply that the initial state affects when the ADM mass 
starts to deviate from its original value. A higher resolution would decrease the discrepancy 
for all initial states and it is expected that all should yield the same thermodynamics. The 
matter Lagrangian tends to zero in the interior in all cases. 
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6.4.2 Extending time evolution in the exterior region: matter Lagrangian and 
exterior RN metric 

The results discussed so far are the robust numerical results of this paper. The exterior region 
has a large outgoing matter wave which can clearly be seen in figure [9] This limits the ability to 
extract unambiguous numerical results in the exterior. Nonetheless, we still want to check two 
things in the exterior, at least qualitatively. We would like to verify that due to the presence 
of the electric field in the exterior, the matter Lagrangian in the exterior does not settle to 
zero in contrast to the Schwarzschild case [U [2] . Secondly, we know the analytical form of the 
RN metric in the exterior. As a consistency check, we would like to verify that numerically 
the metric in the exterior approaches it. 

Unlike the near-horizon region, the exterior does not contain regions with very sharp changes 
in gradients. It is therefore possible to extend the evolution time in the exterior by considering 
only the points outside of a cutoff where one sets some boundary conditions. The cutoff is 
chosen here to be beyond r + at r = 0.5. The outside evolution starts at t = 23.5 and we can 
now relax the time step to At = 10~ 4 . The functions that need to be specified at the cutoff 
are ip, g and the real and imaginary part of x- However, ip and g had already plateaued to a 
constant at the cutoff at t = 23.5; we therefore set them equal to these respective constants on 
the boundary. However, x 1S still evolving at t = 23.5. As time increases, the outgoing wave 
dissipates and we expect the outside to reach the exterior of the RN black hole which contains 
no charge density and a zero stress-energy tensor for the scalar field (there is a non-zero stress- 
energy tensor due to the electric field). This implies that x' should tend to zero outside. Since 
X is zero at infinity, the value of x at the cutoff should tend to zero. All we know is that it 
should approach zero but not how it approaches zero as a function of time. Fortunately, the 
metric is mostly affected by the outgoing wave which is not sensitive to the boundary condition 
at the cutoff. We therefore impose some time-dependent boundary condition for x a t the cutoff 
that leads to a smooth evolution. Figure [11] shows the matter Lagrangian in the exterior at 
very large times. It is important to read this graph properly. The wavy part consists of the 
matter wave which is moving outwards and decreasing in amplitude. What is important is the 
portion before the wave. From t = 40 to t = 120 it remains constant; the matter Lagrangian 
is clearly plateauing towards a non-zero value. This is due to the external electric field. The 
conformal factor ip is plotted in figure [12] together with the analytical form of ip in the exterior 
RN metric. At t = 120, the two curves are very close to each other confirming that the exterior 
is indeed approaching the Reissner-Nordstrom solution. 
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7 Conclusion 

In this paper, we showed that the interior Gibbs free energy during charged collapse stems 
entirely from the gravitational sector and accumulates in the dynamical region inside and near 
the horizon. Numerically, we showed that the interior Lagrangian agrees with the analytical 
expectation Gi n t for the interior Gibbs free energy to within roughly 10% depending on the 
initial state. We observe the formation of a charged shell just inside the horizon with a constant 
electrostatic potential inside that matches &h to within 12%. The work term SW = &h$Q in 
the first law of black hole thermodynamics can be interpreted as the work required to bring a 
charge 5Q to the charged shell at the horizon. 

A consistent dynamical picture of black hole thermodynamics is now emerging from numerical 
studies of gravitational collapse in isotropic coordinates. First, the Helmholtz free energy 
F = E — TS in the Schwarzschild case [HE] makes an interior contribution of Fi n t = —TS and 
the Gibbs free energy G = E—T S — <3? hQ in the charged case makes also an interior contribution 
Gi n t = —TS. The interior contribution in both cases is equal to — TS even though we are dealing 
with different free energies. Secondly, the interior contribution in both cases stems from the 
gravitational sector. Thirdly, in both cases it accumulates inside near the horizon. Fourthly, 
the region where it accumulates is dynamical/non-static. This strongly suggests the following: 
black hole entropy is gravitational entropy and accumulates in the dynamical interior near the 
horizon. 

It is now important to check if gravitational collapse to a Kerr black hole obeys these features. 
Such a study would be a major numerical undertaking. The number of equations to evolve 
would increase but more importantly, a two dimensional spatial grid would now be required, 
increasing massively computation time. The excision technique used in numerical simulations 
of collapse of a rotating neutron star to a Kerr black hole [271 I2H] would not be useful in our 
case. The interior region where the sharp changes in gradients occur in our simulations is 
precisely the region where we need to extract pertinent numerical results. 

Besides the Kerr black hole, it would be interesting to investigate the thermodynamics of the 
black holes recently obtained during numerical studies of collapse of a k-essence scalar field 
[291 [30] . These fields are exotic in the sense that they have a non-canonical kinetic term |31[ [32] . 
It is therefore worthwhile to investigate whether such black holes obey the thermodynamic 
features discussed above. The equations of motion for k-essence scalar fields coupled to gravity 
follow from a Lagrangian and this implies that a thermodynamic study should be possible. 
In particular, one would like to determine if the matter Lagrangian tends towards zero in the 
interior during collapse as with all previous types of matter investigated to date. 
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Ai 

1.85 
1.85 
1.8 



12 

12 



Mtot 

1.364 
1.669 
1.897 



Q 



lol 



M 



-1.234 0.644 
-1.959 0.638 
-2.191 0.754 



Q 

-0.170 
-0.237 
-0.280 



— Lint 

-0.271 
-0.258 
-0.304 



-0.311 
-0.296 
-0.350 



Table 1: Results for profile I 



Table 2: Results for profile II 



Percentage 

discrepancy 

12.6 

12.8 

13.1 



Ai 


A 2 


A3 


Mtot 


Qtot 


M 


Q 


—Lint 


Gint 


Percentage 
discrepancy 


2 


0.11 


0.016 


1.105 


-0.403 


0.706 


-0.110 


-0.315 


-0.349 


9.7 


2 


0.11 


0.032 


1.241 


-0.842 


0.721 


-0.215 


-0.306 


-0.344 


10.9 


2 


0.11 


0.048 


1.488 


-1.357 


0.754 


-0.307 


-0.295 


-0.345 


14.5 



Ai 


Mtot 


M 


— Lint 


Gint 


Percentage 
discrepancy 


1.50 


1.102 


0.334 


-0.295 


-0.334 


11.5 


1.55 


0.974 


0.279 


-0.244 


-0.279 


12.5 


1.60 


0.867 


0.236 


-0.203 


-0.236 


14.0 



Table 3: Results for profile III 



Number 


of 


spatial 


At 


steps 








1250 






0.0008 


2500 






0.0004 


5000 






0.0002 


10000 






0.0001 


20000 






0.00005 



Time of deviation Time of deviation 
of 1% of Mt t of 5% of M^t 



14.3 
16.4 
18.3 
20.1 
21.7 



17.9 
19.4 
20.9 
22.3 
23.7 



Table 4: Evolution time as a function of gridsize for profile II with Ai 
A 3 = 0.032 
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Figure 1: The position of outgoing radial light-like geodesies is shown as a function of time. 
They are released at t = 12. The radius r + where N = is also shown as a function of time 
(the top black line). The null geodesies are unable to cross r + (see also [23]). 
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Figure 2: Space profile of the metric field ip for different times. 
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Figure 3: Space profile of the lapse function N for different times. r + is the radius where N 
crosses zero at late times. TV is negative in the interior and approaches zero from below. 
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Figure 4: Space profile of the absolute value of the matter field x f° r different times. Note the 
propagation of an outgoing matter wave in the exterior region. 
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Figure 5: Potential a for different times. Inside the horizon, a is almost constant, consistent 
with a charged shell configuration at the apparent horizon. The change at large r in the profile 
of a corresponds to the position of the outgoing wave. 
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Figure 6: The charge density at different times. As the system evolves, it becomes more and 
more concentrated towards the inside of the apparent horizon. At late times one has a charged 
shell in the vicinity of the horizon. 
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Figure 7: The gravitational Lagrangian accumulation at different times. Note the negative 
dip in the thin shell just inside and near the apparent horizon. There is a disturbance in the 
outgoing wave region. 
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Figure 8: The time-derivative of the metric function if), ip is static in the exterior and its 
time-derivative is decreasing towards zero in the interior volume. However, it is not static 
inside and near the horizon, precisely the region where the interior free energy accumulates. 
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Figure 9: The matter Lagrangian accumulation at different times. Note that there is no 
disturbance at the apparent horizon at late times. The disturbance is around the outgoing 
wave. 
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Figure 10: The minimum in the Lagrangian accumulation and the interior Gibbs free energy 
Gi n t as a function of time. The interior Gibbs free energy is given by — r + . 
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Figure 11: The matter Lagrangian accumulation in the exterior plotted at different times. 
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Figure 12: The conformal factor ijj at t = 120 is compared to its theoretical expectation in the 
exterior. 
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